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A  novel  modeling  scheme  of  SOFC  anode  based  on  the  stochastic  reconstruction  technique  and  the  Lat¬ 
tice  Boltzmann  Method  (LBM)  is  proposed  and  applied  to  the  performance  assessment  and  also  to  the 
optimization  of  anode  microstructures.  A  cross-sectional  microscopy  image  is  exploited  to  obtain  a  two- 
dimensional  phase  map  (i.e.,  Ni,  YSZ  and  pore),  of  which  two-point  correlation  functions  are  used  to 
reconstruct  a  three-dimensional  model  microstructure.  Then,  the  polarization  resistance  of  the  recon¬ 
structed  anode  is  obtained  by  the  LBM  simulation.  The  predicted  anodic  polarization  resistance  for  a 
given  microstructure  and  its  sintering  temperature  dependence  are  in  good  agreement  with  the  literature 
data.  Three-dimensional  distributions  of  potential  and  current  can  be  obtained,  while  and  the  effect  of 
working  temperature  is  discussed.  The  proposed  method  is  considered  as  a  promising  tool  for  designing 
SOFC  anodes. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  considered  to  be  one  of  the 
most  powerful  candidates  for  the  future  energy  systems,  because 
of  their  high  efficiency  and  fuel  flexibility  [1  ].  In  addition,  their  high 
temperature  waste  heat  can  be  utilized  in  a  heat  engine,  achieving 
an  overall  efficiency  over  60%  [2],  and  also  for  fuel  reforming.  On 
the  other  hand,  improved  performance  at  lower  temperature  oper¬ 
ation  is  desirable  for  cost  reduction  and  system  durability.  Recently, 
anode-supported  SOFCs  with  thin  electrolyte  have  attracted  inter¬ 
est  [3],  because  their  low  ohmic  loss  enables  lower  temperature 
operation. 

Generally,  porous  cermets  (such  as  Ni-YSZ  composite)  are  used 
as  anode  materials  in  order  to  match  the  thermal  expansion  coef¬ 
ficient  with  that  of  electrolyte  and  to  extend  effective  reaction 
sites  (three-phase  boundary,  TPB).  Anode  microstructure  parame¬ 
ters  (porosity,  grain  diameter,  etc.)  are  known  to  have  significant 
effects  upon  the  cell  performance  and  durability  [4,5].  For  the 
detailed  evaluations  of  these  parameters,  numerical  modeling  of 
anode  microstructures  has  been  strongly  demanded  [6-10]. 

Conventional  electrode  modeling  is  mainly  classified  into  three 
types,  i.e.,  (i)  thin-film  model  (or  pore  model)  [6],  (ii)  random  resis¬ 
tor  network  model  [7,8],  and  (iii)  random  packing  sphere  model 
[9,10].  In  the  first  model,  the  network  of  ionic  conductor  is  treated 
as  tree-like  protrusion  covered  with  electronic  conductors.  In  the 
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second  model,  particles  of  ionic  and  electronic  conductors  are  ran¬ 
domly  distributed  by  a  Monte  Carlo  technique,  and  modeled  as  a 
resistor  network.  Then,  the  overpotential  in  modeled  electrode  is 
obtained  by  solving  Kirchoff  s  law  in  this  resistor  network.  The  third 
model  is  an  advanced  model  of  classical  electrode  models  [  11  ] :  ionic 
and  electronic  conductivities  and  TPB  density  are  calculated  by 
applying  the  percolation  theory  [12]  and  the  coordination  number 
theory  [13]  of  randomly  packed  spheres.  Most  of  these  models  rely 
on  the  macroscopic  homogeneity  of  composite  electrode  media, 
and  their  consistency  with  the  actual  electrode  structures  experi¬ 
mentally  observed  [14]  has  not  been  assured. 

Recently,  detailed  microscopic  (grain-scale)  simulation  of  SOFC 
electrodes  has  been  paid  more  attention  [15,16].  Joshi  [15]  per¬ 
formed  detailed  multi-component  diffusion  simulation,  while 
Abbaspour  et  al.  [16]  solved  the  species-transport  problem 
coupled  with  electrochemical  reaction  in  modeled  electrodes. 
However,  their  electrode  model  is  only  two-dimensional  being 
over-simplified  as  somewhat  artificial  configuration.  Recently, 
Wilson  et  al.  [17]  performed  the  detailed  measurement  of  3D 
microstructure  of  anodes  using  FIB-SEM,  and  solved  the  mass 
transport  in  the  measured  microstructure  by  FEM.  However  it 
requires  huge  efforts  for  measurement  and  data  processing,  and 
the  effect  of  microstructures  upon  electrochemical  activity  remains 
unclear. 

Mukherjee  and  Wang  [  18  ]  proposed  a  detailed  modeling  scheme 
of  species  and  charge  transport  in  a  polymer  electrolyte  fuel  cell 
(PEFC)  catalyst  layer.  The  realistic  porous  microstructure  of  the 
catalyst  layer  was  reconstructed  by  using  the  stochastic  reconstruc¬ 
tion  scheme  originally  proposed  by  Quiblier  [19].  This  modeling 
scheme  could  offer  better  understanding  of  PEFC  microstructures 
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Nomenclature 

C  molar  density  (mol  m3) 

D  diffusion  coefficient  (m2  s-1 ) 

E  cost  function  of  reconstruction 

fi  velocity  distribution  function 

fteq  equilibrium  velocity  distribution  function 

i  current  density  (A  m-2 ) 

z'o  exchange  current  density  per  TPB  length  (A  m-1 ) 

z'reac  reaction  current  density  (A  nrr3 ) 

L  electrode  thickness  (pun) 

Lj(r)  lineal-path  function 

/j(r)  number  of  line  elements 

Ltpb  TPB  density  per  unit  volume  (m-2 ) 

Npixei  number  of  pixels 

p  pressure  (Pa) 

Rij(r)  two-point  correlation  function 

t*  relaxation  time 

wz  reaction  production 

Z(x)  phase  function 

Greek  symbols 
s  volume  fraction 

0  potential  (V) 

r  Thiele  modulus 

p  over  voltage  (V) 

p  resistivity  (£2  m) 

a  conductivity  (S  nrr 1 ) 

r  tortuosity  factor 

Subscripts  and  superscripts 

eff  effective  value  in  porous  media 

elec  electron 

ion  ion 

H2  hydrogen 

H20  vapor 


and  their  effects  upon  cell  performance.  However,  their  catalyst 
layer  model  assumed  only  two  phases,  z.e.,  gas  phase  and  mixed 
electrolyte/electrode  phase,  and  empirical  Bruggman  correlation 
was  employed  for  the  solid  phase  conduction.  In  SOFC  anodes, 
the  ionic  conductor  (YSZ),  the  electronic  conductor  (Ni)  and  the 
gas  phase  conductor  (pore)  are  comparable  in  size  and  volume 
fraction,  so  that  such  two-phase  assumption  cannot  be  readily 
verified. 

In  this  study,  anode  microstructures  are  reconstructed  by  using 
the  reconstruction  scheme  for  multiphase  media  proposed  by 
Yeong  and  Torquato  [20].  Then,  the  anodic  overpotential  of  the 
reconstructed  structures  is  obtained  by  using  the  Lattice  Boltz¬ 
mann  Method  (LBM)  [21].  The  effects  of  sintering  temperature 
upon  anode  microstructures  and  cell  performances  are  evaluated. 
Furthermore,  the  effects  of  working  temperature  on  activation 
polarizations  as  well  as  on  local  three-dimensional  distributions 
of  potential  and  current  are  discussed. 

2.  Sample  preparation 

Composite  anode  samples  were  prepared  by  mixing  8mol% 
YSZ  and  nickel  oxide  (NiO)  powders  (AGC  Seimi  Chemical  Co. 
Ltd.,  Japan,  averaged  grain  diameter  =  1.5  p.m).  Mixing  ratio  of  NiO 
and  YSZ  was  60:40  wt%,  and  the  Ni-YSZ  volume  ratio  resulted  as 
42.4:57.6  after  reduction.  Samples  were  molded  into  rods  by  extru¬ 
sion  molding  (Kankyo  Ceramics  Research  Co.  Ltd.,  Japan).  Three  sets 


Fig.  1.  Height  image  of  a  1400°C-sintered  sample  measured  by  a  scanning  probe 
microscope. 


of  samples  were  pre-sintered  at  1000  °C  for  3h,  and  sintered  for 
3h.  Sintering  temperature  was  changed  as  1300,  1350  or  1400  °C 
for  evaluating  the  effects  of  sintering  temperatures.  Then,  sam¬ 
ple  anodes  were  reduced  in  67%  H2  diluted  with  N2  at  750 °C  for 
10  h.  The  ramp  rate  was  10  °C  min-1  in  all  heating  and  cooling  pro¬ 
cesses.  A  porosity  of  a  1200°C-sintered  sample  was  measured  by 
a  mercury  porosimeter  (Autopore  IV9500,  Shimadzu,  Japan).  The 
porosities  of  1300, 1350  and  1400  °C-sintered  samples  extrapolated 
from  their  shrinkage  ratios  were  0.450,  0.393  and  0.335,  respec¬ 
tively. 

After  reduction,  samples  were  infiltrated  with  epoxy  to  avoid 
smearing  of  soft  Ni,  and  polished  with  a  diamond  paste.  For  the 
stochastic  reconstruction  described  later,  highly  flat  cross-sectional 
surfaces  should  be  obtained.  In  order  to  evaluate  the  surface 
roughness,  cross-section  of  a  polished  sample  is  measured  with  a 
scanning  probe  microscope  (SPM9600,  Shimadzu,  Japan;  512  x  512 
pixels,  dynamic  mode  AFM).  As  shown  in  Fig.  1,  the  surface  rough¬ 
ness  was  less  than  80  nm,  and  considered  to  be  negligibly  small 
compared  to  the  grain  size. 

Cross-sectional  optical  images  of  polished  samples  were 
obtained  with  a  confocal  laser  microscope  (OLS3000,  Shimadzu, 
Japan;  1024  pixel  x  768  pixel,  1  pixel  =  41  nm).  The  original 
microscopy  images  are  shown  in  Fig.  2.  Ni  (white),  YSZ  (gray) 
and  pore  (black)  can  be  clearly  distinguished.  Grain  growth 
under  high-temperature  sintering  was  observed,  especially  for 
Ni. 


3.  Stochastic  reconstruction 

3.1  2D  image  processing 

Three  phases  in  raw  images  (Fig.  2)  can  be  easily  distinguished 
by  their  brightness  values  (see  Fig.  3).  To  sharpen  grain  boundaries, 
a  weighted  differential  filter  was  applied.  Then,  each  phase  was 
distinguished  by  applying  the  scheme  of  Lee  et  al.  [22].  Threshold 
of  brightness  value  was  determined  so  that  the  volume  fractions 
in  the  filtered  image  match  with  the  values  calculated  from  the 
porosity  and  the  Ni-YSZ  volume  ratio.  A  phase-distinguished  image 
of  the  1400  °C-sintered  sample  is  shown  in  Fig.  4,  where  clear  phase 
differentiation  is  observed. 
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Fig.  2.  Original  cross-sectional  images  of  sample  anodes  sintered  at  (a)  1300  °C,  (b)  1350  °C  and  (c)  1400  °C. 


Fig.  3.  Brightness  histogram  of  a  1400°C-sintered  sample. 


3.2.  2D  image  analysis 

From  the  phase-distinguished  2D  images,  a  phase  function,  Z(x), 
is  defined  at  each  pixel  (x)  as  follows: 

{0  (forx  g  pore,  black) 

1  (for xeYSZ,  gray)  .  (1) 

2  (forx g  Ni,  white) 


ficient  imaging  areas,  but  the  error  was  relatively  small  (less  than 
5%).  As  shown  in  Fig.  5,  the  YSZ-YSZ  correlation  seems  to  be  more 
uniform  than  that  of  Ni-Ni,  and  this  represents  that  YSZ  grains  are 
better  dispersed  in  cermet  anodes. 


3.3.  Stochastic  reconstruction 

Stochastic  reconstruction  was  originally  proposed  byjoshi  [23], 
and  later  extended  into  three-dimensions  by  Quiblier  [19].  How¬ 
ever,  their  schemes  based  on  the  Gaussian  filtering  are  relatively 
poor  in  accuracy,  and  also  not  suitable  for  multiphase  media.  There¬ 
fore,  the  iterative  scheme  of  Yeong  and  Torquato  [20]  is  applied  in 
this  study.  The  grid  length  Ax  is  0.178  \xm  in  each  direction,  which 
is  around  1/10  of  the  averaged  grain  diameter.  This  grid  size  is  con¬ 
sidered  to  be  small  enough  from  our  preliminary  calculation.  The 
whole  computational  domain  size  is  150  grids  (26.7  p,m  in  space) 
cubic,  and  considered  to  be  thick  enough  compared  to  the  effective 
anode  thickness  (~20  |jim)  [24].  Periodic  boundary  conditions  are 
imposed  in  each  direction.  At  first,  each  voxel  is  assigned  randomly 
with  three  phases  according  to  their  volume  fractions.  Then,  recon¬ 
struction  is  carried  out  to  minimize  the  cost  function  E  defined  as 
follows: 


Using  this  phase  function,  the  two-point  correlation  (or  probability) 
function  Ry  between  phases  i  and  j  (ij  =  0-2)  is  defined  as  follows: 


%(<-)= 


3(Z(x),i)3(Z(x  +  r),j) 

«(Z(x),  T) 


(2) 


where  8{i}j)  represents  the  Kronecker’s  delta  function.  Fig.  5  shows 
the  two-point  correlation  of  Ni-Ni  and  YSZ-YSZ.  These  correlation 
functions  did  not  fully  converge  on  volume  fractions  due  to  insuf- 


Fig.4.  Phase-distinguished  image  of  a  1400  °C-sintered  sample  (Ni:  white,  YSZ:  gray, 
pore:  black). 


{Rf(r)-Rf(r)}2dr, 


(3) 


where  R?P  and  RlP  are  the  reconstructed  and  the  original  2D  two- 
point  correlation  functions,  respectively.  Isotropic  porous  media  are 
assumed,  and  the  integral  dr  is  performed  only  in  the  directions  of 
the  principal  axes  in  order  to  reduce  the  computational  load.  The 
interval  of  integration  r0  was  set  as  7  p,m. 

Starting  from  a  random  configuration,  two  voxels  of  different 
phases  are  exchanged  at  each  iteration  step  only  when  the  cost 
function  E  decreases.  The  reconstruction  is  terminated  after  the  cost 
function  remained  constant  for  50  Monte  Carlo  steps.  Then,  resulted 
porous  structures  are  modified  by  using  the  surface  curvature 
modification  scheme  [25]  in  order  to  avoid  physically  impossible 
structures  such  as  “solid  phase  floating  in  a  pore”. 


3.4.  Validity  of  reconstructed  structure 

Reconstructed  3D  microstructures  of  three  samples  are  shown 
in  Fig.  6.  Note  that,  in  the  figures,  small  parts  of  reconstructed  struc¬ 
tures  (100  grid  cubic)  are  smoothed  out  for  eye-friendliness.  Good 
agreement  is  achieved  between  reconstructed  and  original  two- 
point  correlation  functions  as  shown  in  Fig.  7(a),  with  the  residuals 
of  Eq.  (3)  being  less  than  1.0%.  In  order  to  check  the  validity  of 
the  reconstructed  structures,  the  lineal-path  function  If(r)  of  the 
reconstructed  structure  is  compared  with  that  of  the  original  2D 
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Fig.  6.  Reconstructed  anode  microstructures  sintered  at  (a)  1300  °C,  (b)  1350  °C  and  (c)  1400  °C.  Green:  Ni,  blue:  YSZ,  transparent  gray:  pore.  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article). 


Fig.  7.  (a)  Two-point  correlation  functions  and  (b)  lineal-path  functions  of  original  2D  image  and  reconstructed  3D  structure  (1400°C-sintered  sample). 


images. 


C(r)  = 


^pixel 


(4) 


In  Eq.  (4),  Z2(r)  represents  the  number  of  line  segments  of  length  r 
lying  in  phase  i,  while  Npixel  represents  the  number  of  total  pixels. 
As  shown  in  Fig.  7(b),  reconstructed  lineal-path  functions  match 
well  with  those  of  the  original  2D  image. 


4.  Electrode  performance  simulation  scheme 


ing  condition,  current  collector  (CC)  and  YSZ  electrolyte  (EL)  layers 
are  added  to  the  boundaries  at  x  =  0m  and  x  =  26.7  pm  (150 Ax), 
respectively.  The  thicknesses  of  the  CC  and  the  EL  layers  are  1.78  pm 
(lOAx)  and  3.56  pm  (20Ax),  respectively.  For  simplicity,  binary 
gas  composition  of  H2/H20  is  assumed,  and  equimolar  diffusion 
assumption  [26]  is  employed, 


D  = 


-°»2  +  ~jfDH20- 


(5) 


4.1  Simulation  model 


The  diffusion  coefficient  of  each  species  i  is  determined  as: 


The  anodic  overpotential  of  reconstructed  anodes  is  calculated 
from  detailed  numerical  simulations.  A  schematic  of  computational 
domain  is  shown  in  Fig.  8.  In  order  to  model  the  actual  anode  operat- 


Di  = 


Um  +D,J 


(6) 
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Constant  ionic  current  density 


Electrolyte  layer 
anode  layer 
^urrem  collector  layer 
Fig.  8.  Schematic  of  a  computational  domain. 


Constant  gas 
&  electronic  current  density 


In  Eq.  (6),  Dim  and  Dik  represent  binary  and  Knudsen  diffusion  coef¬ 
ficients,  respectively,  being  given  as: 


Pore 


RT  . 

*  2FPi,eac  * 


where  r\  is  the  overpotential  </>eiec  -  0ion.  The  lineal  exchange  cur¬ 
rent  density  i0  is  defined  from  the  porous  Ni  anode  experiments  of 
De  Boer  [30]  as: 


Dim  =  5.956 


1/2 


T3/ 2 

D  ’ 


(7) 


(8) 


„  2/8  RT\V2 

Dik-3{m;)  rp- 

Gas  properties  are  summarized  in  Table  1,  and  the  collision  integral 
is  represented  as: 

-0.1814 


£2D  =  1.1336 


(?)' 


(9) 


The  rp  in  Eq.  (8)  represents  the  mean  pore  radius,  and  was  set  as 
0.75  |jim  for  all  simulated  cases. 

Temperature  is  considered  to  be  uniform,  and  uniform  current 
is  assumed  at  x  =  -1.78  pan  (-10Ax)  and  x  =  30.26  jxm  (170 Ax).  The 
material  properties  [27,28]  are  summarized  in  Table  1. 

4.2.  Governing  equations 


The  governing  equations  in  the  pore,  Ni  and  YSZ  phases  are  given 
as  the  diffusion  equation  of  hydrogen  gas  and  the  charge  transfer 
equations  for  electron  and  ion,  respectively.  They  are  written  as 
follows: 


V(DVC)  =  — -rrlr, 

2  F 

V(cre  lec^^elec)  —  —  heac, 

^(^ion^^ion)  —  heac- 


(10) 
(11) 
(12) 

The  RHSs  in  Eqs.  (10)-(12)  are  the  reaction  production  terms, 
and  the  local  reaction  current  ireac  is  calculated  by  using  the 
Butler-Volmer  equation  [29]: 


:  =  io^TPB 


-N 

/  F  \] 

rpiRf"j 

“eXP  {-RfV\ 

(13) 


Table  1 

Material  properties  used  in  simulation 


Properties 

Value  or  expression 

Electronic  conductivity  crelec  (S  m_1 )  [27] 

Ionic  conductivity  cri0n  (S  m-1 )  [27] 

3.37  x  106  -  1065.37 

3.37  x  104  exp(- 10300/7) 

Substance 

M  (10-3  kg  moP1) 

£(A)  8//<d<) 

h2 

h2o 

2.016 

18.015 

2.93  37 

2.65  356 

io  =  31.4  X  ph°'0Vh24o°  exp  (-1'83;104)  .  (14) 

The  three-phase  boundary  length  per  unit  volume,  LTPB,  in  Eq.  (13) 
is  assumed  to  be  20%  smaller  than  the  value  directly  calculated  from 
the  cubic  voxel  perimeter  of  the  reconstructed  structure.  Geomet¬ 
rically,  the  actual  TPB  length  should  lie  between  57.7%  and  100%  of 
the  cubic  voxel  perimeter.  In  the  present  study,  the  voxel  perimeter 
was  simply  reduced  by  20%  as  a  representative. 

4.3.  Numerical  scheme 

The  LBM  is  used  to  solve  Eqs.  (10)-(12)  in  each  phase.  The  LB 
equation  with  the  LBGK  model  in  the  collision  term  is  written  as 
follows: 

/i'(x  +  qAt,  t  +  At)  =fi(x,  t)  -  t )  -f-q(x,  t)]  +  w,At.  (15) 

In  Eq.  (15),  represents  the  velocity  distribution  function  of  gas, 
electron  or  ion  with  velocity  cz  in  the  ith  direction,  and  is  the 
Maxwellian  local  equilibrium  distribution.  For  the  3D  LBM  simula¬ 
tion,  D3Q15  (i  =  1 — 15)  or  D3Q19  (i  =  1—19)  models  are  commonly 
used.  However,  it  has  been  shown  that,  in  case  of  simple  diffu¬ 
sion  simulation,  D3Q6  (i  =  1  —6)  model  can  be  used  with  a  slight 
loss  of  accuracy  [31  ].  So  the  D3Q6  model  is  used  in  this  work.  The 
relaxation  time  t*  is  0.99  and  fixed  for  all  simulations. 

Periodic  boundary  conditions  are  assumed  in  the  spanwise  (y 
and  z)  directions.  At  the  CC  surface,  constant  gas  composition  is 
assumed.  Constant  electronic  and  ionic  current  flux  conditions  are 
imposed  on  the  CC  and  EL  boundaries.  A  no-flux  boundary  condi¬ 
tion  is  imposed  on  the  boundary  of  each  phase  in  the  porous  media 
by  applying  the  halfway  bounceback  scheme  with  second-order 
accuracy  [32]. 

The  production  term  w*  in  the  RHS  is  calculated  from  Eqs. 
(10)-(14)  using  a  scheme  similar  to  that  of  Mukherjee  and  Wang 
[18].  A  schematic  of  reaction  calculation  is  shown  in  Fig.  9.  The 
reaction  current  at  TPB  (t  +  l/2,j  + 1/2,  k)  is  calculated  as: 

ireac  (^i  +  ^  ^ »  k'j  =  io  x  f'TPB  ^  ,j +  l<^j 

{exp((2 F/RT)q)  -  exp (-(F/RT)^)} 

X  Ax3 

(16) 

4  —  0elec (bj\  —  0ion(i  +  1  ^)-  (17) 
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Fig.  10.  (a)  Conductivity  of  solid  phases  and  (b)  tortuosity  factor  of  gas  phase. 


Table  2 

Simulation  conditions  in  Section  5.1 


Properties 

Value  or  expression 

Operating  temperature  T  (°C) 

1000 

Pressure  p  (Pa) 

1.013  x  105 

Fuel  composition  (H2:H20)  (mol%) 

97:3 

Current  density  i  (A cm-2) 

■sr 

O 

X 

o 

5.  Computational  results 

5.1.  Electrode  performance 


The  calculated  electronic  and  ionic  conductivities,  crfJL  and  of£, 

’  elec  ion’ 

and  gas  phase  tourtuosity  factor,  r,  are  shown  in  Fig.  10.  With 
increasing  the  sintering  temperature,  the  solid  phase  conductivities 
are  improved  while  gas  phase  connectivity  is  reduced.  The  calcu¬ 
lated  gas  phase  tortuosity  factors  are  in  good  agreement  with  the 
literature  data  (r  =  2.5-4.0)  [33]. 

The  reconstructed  anode  performances  are  compared  with  the 
experimental  data  [34].  The  simulation  conditions  are  summarized 
in  Table  2.Since  the  experimental  data  was  taken  from  the  electro¬ 
chemical  impedance  spectroscopy  under  zero  bias  current,  only  a 
small  value  of  current  i  =  lx  10-4  A  cm-2  was  assumed  in  this  case. 
The  temperature  dependence  of  the  anodic  polarization  resistance 
is  shown  in  Fig.  11.  It  is  found  that  the  best  anode  performance 
is  achieved  by  the  sample  sintered  at  1400  °C.  This  tendency  is 
in  good  agreement  with  the  experimental  results  in  literatures 
[34,35],  especially  with  the  results  of  Primdahl  et  al.  [34],  whose 
sample  condition  is  similar  to  the  present  case. 

This  sintering  temperature  dependence  can  be  explained  by 
Thiele  modulus  r  (effectiveness  criterion)  [36],  which  can  be 


Fig.  12.  Sintering  temperature  dependence  of  the  Thiele  modulus  r. 


obtained  by  linearizing  and  coupling  Eqs.  (11 )— (13): 


2  =  3»qLtpbF«  +  <c) 
RT 


(18) 


Thiele  modulus  represents  the  ratio  of  ohmic  resistance  to  the 
activation  overpotential.  Generally,  ohmic  resistance  becomes 
dominant  and  reactive  electrode  thickness  becomes  inversely  pro¬ 
portional  to  r  when  r>  3  [36].  For  such  cases,  the  order  of  the 
anode  polarization  loss  can  be  expressed  as: 


Tj  oc 


I'TPB  *0 


(19) 


In  Eq.  (19),  electronic  resistivity  is  neglected  because  it  is  much 
smaller  than  the  ionic  resistivity. 

In  Fig.  12,  the  sintering  temperature  dependence  of  Thiele  mod¬ 
ulus  r  is  shown.  Thiele  modulus  r  of  each  sample  is  larger  than 
3,  and  thus  ohmic  (ionic)  resistance  is  dominant.  Eq.  (19)  reveals 
that  the  anode  polarization  loss  is  dependent  on  the  ionic  resis¬ 
tivity  pf*fn  and  the  TPB  density  Ltpb-  The  calculated  ionic  resistivity 


Fig.  11.  Effect  of  sintering  temperature  upon  the  anodic  polarization  resistance 
(comparison  with  literature  data). 


Fig.  13.  Sintering  temperature  dependence  of  the  TPB  density  and  the  ionic  resis¬ 
tance. 
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Fig.  14.  (a)  Phase,  (b)  potential  and  (c)  current  distribution  at  anode-electrolyte  interface  of  a  1400  °C-sintered  sample.  Phases  in  (a)  are;  white:  Ni,  light  gray:  YSZ,  dark  gray: 
pore,  black:  disconnected  phases,  red  lines:  TPB.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


10  Jim 


B 40000 
—  ([A/m2] 


Fig.  15.  Temperature  dependence  of  ionic  current  distribution  at  anode-electrolyte  interface;  (a)  600  °C,  (b)  800  °C  and  (c)  1000  °C. 


Table  3 

Simulation  conditions  in  Section  5.2 


and  the  TPB  densities  (Ltpb)  are  shown  in  Fig.  13.  With  increas¬ 
ing  the  sintering  temperature,  the  ionic  resistivity  decreases  while 
the  TPB  density  increases,  which  leads  to  smaller  polarization 
loss. 

Also,  the  TPB  densities  in  the  modeled  anodes  are  found 
to  be  around  2.4-2.8  |jinrr2.  This  predicted  TPB  densities  are 
somewhat  smaller  than  the  value  obtained  from  FIB-SEM 
measurement  4.28  jim-2  [37].  Flowever,  the  order  of  magni¬ 
tude  is  the  same  although  they  are  obtained  from  different 
samples. 

5.2.  Local  potential  and  current  distributions 

The  effects  of  working  temperature  on  local  distributions  of 
potential  and  current  are  investigated.  Simulation  conditions  are 
summarized  in  Table  3.  In  order  to  investigate  a  high  power  den¬ 
sity  operating  condition,  the  current  density  is  set  as  i  =  1 A  cm-2. 
In  Fig.  14,  the  phase,  current  and  potential  distributions  in 
an  anode-electrolyte  interface  at  x  =  26.7|jim  (150 Ax)  of  the 
1400°C-sintered  sample  are  shown.  The  y-z  plane  is  in  con¬ 
tact  with  the  electrolyte  (YSZ)  surface,  and  thus  the  pore-Ni 
interfaces  (red  lines)  represent  TPBs.  As  shown  in  Fig.  14,  the 
potential  and  current  distributions  are  considerably  inhomoge¬ 
neous. 

The  predicted  current  distributions  for  different  working  tem¬ 
peratures  are  shown  in  Fig.  15.  With  decreasing  the  working 
temperature,  the  current  density  around  TPBs  is  significantly 
increased  (see  Fig.  14(a)).  In  Fig.  16,  the  ionic  current  density 
distribution  in  the  x  direction  is  shown.  Under  600  °C  opera- 


Fig.  16.  Effect  of  working  temperature  upon  the  effective  anode  thickness. 


tion,  the  reactive  electrode  thickness  becomes  very  thin  (<5  p,m), 
and  this  fact  means  that  the  reaction  is  confined  in  this  thin 
region  due  to  relatively  large  ionic  resistivity.  Thus,  the  ionic 
current  flows  mainly  in  the  vicinity  of  TPB.  Even  though  Thiele 
modulus  r  decreases  with  the  operating  temperature  as  shown 
in  Fig.  17,  non-linear  dependence  of  activation  overpotential  on 


Fig.  17.  Working  temperature  dependence  of  the  anodic  overpotential  and  Thiele 
modulus  r. 
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current  makes  the  ionic  resistance  dominant  for  this  low  tem¬ 
perature  and  large  current  density  condition.  In  other  words,  the 
linearized  Butler-Volmer  equation  in  Thiele  modulus  r  becomes 
invalid  because  exchange  current  density  i0  becomes  small  for 
such  conditions.  In  the  case  of  600  °C  operation,  the  activa¬ 
tion  resistance  calculated  from  Eq.  (13)  is  much  smaller  than 
the  ionic  resistivity,  making  the  reactive  anode  thickness  quite 
thin. 

6.  Conclusions 

A  novel  modeling  scheme  of  SOFC  anodes  based  on  the  stochas¬ 
tic  3D  structure  reconstruction  scheme  and  the  LBM  simulation  is 
proposed  and  evaluated  for  the  prediction  of  the  effects  of  anode 
microstructures  on  the  total  anode  performance.  Cross-sectional 
2D  microscopy  images  of  anodes  are  processed  to  obtain  phase- 
distinguished  images,  and  3D  anode  microstructure  models  are 
reconstructed  from  the  two-point  correlation  functions  of  the  2D 
images.  Good  reconstruction  is  achieved  by  the  present  scheme, 
and  residuals  in  two-point  correlation  functions  and  lineal-path 
functions  are  negligibly  small. 

The  anodic  overpotential  of  reconstructed  anodes  is  calculated 
by  coupling  the  mass  and  charge  transfer  and  the  electrochemi¬ 
cal  reaction  using  LBM,  and  the  effect  of  anode  microstructures  on 
the  cell  performance  is  discussed.  It  is  found  that  the  optimal  sin¬ 
tering  temperature  should  be  around  1400  °C,  where  a  good  ionic 
network  and  electrochemical  activity  are  obtained.  This  result  is  in 
good  agreement  with  the  existing  literature  data. 

By  using  this  scheme,  three-dimensional  distributions  of  current 
and  potential  are  obtained,  and  the  effects  of  working  temperature 
on  local  potential  and  current  distributions  are  discussed.  With  the 
decrease  of  working  temperature,  the  current  concentrates  near 
TPB  of  anode/electrolyte  interface  and  the  reactive  anode  thickness 
becomes  thinner  (<5  p,m).  This  is  attributed  to  the  dominance  of 
ionic  resistance  over  activation  overpotential. 

This  scheme  gives  light  to  better  understanding  of  microscale 
phenomena  and  the  optimal  design  of  microstructures  in  SOFC 
anodes. 

The  present  scheme  requires  validation  with  the  directly  mea¬ 
sured  3D  data,  such  as  FIB-SEM  measurement  in  Ref.  [17].  And  this 
remains  as  a  future  challenge. 
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